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1.  Introduction 

Suppose  one  injects  a  drug  into  an  animal.  What  is  the  time  course  of  the  drug 
concentration  in  the  blood  and  in  the  various  organs  of  the  body?  Although  pre¬ 
liminary  attempts  to  set  up  mathematical  models  of  this  problem  were  made  by 
Teorell  [1],  [2]  and  models  of  certain  special  cases  have  been  considered  by  von 
Schrotter  [3],  Smith  and  Morales  [4],  and  Morales  and  Smith  [5],  [6],  it  con¬ 
tinues  to  challenge  the  biologist  and  mathematician.  Much  of  the  recent  theo¬ 
retical  work  in  this  area  has  been  concerned  with  the  analysis  of  radioactive 
tracer  experiments  on  compartmentalized  systems  assumed  to  be  in  a  steady 
state  [7],  [8].  This  leads  to  the  usual  sets  of  simultaneous  linear  differential 
equations.  We  wish  to  consider  the  more  complex  kinetic  problem.  In  particular, 
we  wish  to  consider  the  kinetics  of  distribution  of  an  injected  compound  assum¬ 
ing  that  it  exchanges  between  capillary  plasma  and  extracellular  space  by  dif¬ 
fusion,  enters  the  intracellular  space  by  diffusion  and  perhaps  by  some  active 
transport  process  and  can  react  with  some  constituent  (say  an  enzyme)  of  the 
intracellular  space.  Although  the  motivation  for  studying  this  problem  was  the 
desire  to  obtain  a  theoretical  framework  for  the  analysis  of  results  with  various 
agents  used  in  cancer  chemotherapy,  the  problem  is  of  more  general  importance 
for  physiology  and  pharmacology. 

2.  Anatomical  and  physiological  considerations 

Any  mathematical  model  of  such  a  complex  process  must  be  firmly  imbedded 
in  the  known  anatomical  and  physiological  background  if  it  is  to  approach 
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reality.  Let  us  therefore  review  briefly  the  basic  anatomy  and  physiology  of  the 
capillary  circulation. 

A  tissue  or  region  of  an  organ  which  we  will  be  concerned  with  consists  of 
many  living  cells.  Most  of  the  cells  of  any  particular  tissue  have  roughly  similar 
shapes  and  sizes.  Between  the  cells  is  a  watery  fluid,  the  extracellular  fluid. 
Each  cell  is  bounded  by  a  cell  membrane  which  is  lipoprotein  in  nature.  The 
lipid  layer  of  the  cell  membrane  acts  as  a  barrier  to  the  diffusion  of  polar  mole¬ 
cules  into  the  cell:  so  much  so  that,  for  polar  molecules,  diffusion  across  the  cell 
membrane  is  orders  of  magnitudes  slower  than  diffusion  across  an  aqueous  layer 
of  comparable  thickness.  The  blood  supply  to  such  a  region  enters  via  arteries 
which  branch  and  eventually  give  rise  to  the  precapillary  arterioles.  The  arteri¬ 
oles  give  rise  to  many  small  tubules  whose  walls  consist  of  a  single  thin  layer  of 
flat  cells.  These  are  the  capillaries  which  pass  through  the  extracellular  space 
and  finally  join  to  form  venules  which  combine  to  form  the  veins  which  carry 
the  blood  away  from  the  region.  According  to  Zweifach  [9],  [10],  there  are  pre¬ 
ferred  channels,  A-V  bridges,  connecting  arterioles  and  venules  and  many  of  the 
capillaries  arise  from  these  A-V  bridges.  These  A-V  bridges  have  muscle  cells  at 
various  points  of  their  walls  which  presumably  function  in  the  regulation  of  the 
capillary  circulation.  With  low  blood  flows  to  a  region,  the  main  flow  is  through 
the  A-V  bridge.  As  the  flow  increases,  more  and  more  of  the  flow  passes  through 
the  capillaries  arising  from  the  A-V  bridge.  Thus  increased  blood  flow  is  to  a 
great  extent  accommodated  by  an  opening  of  more  channels  and  not  simply  by 
an  increased  flow  rate  in  a  fixed  number  of  channels.  It  follows  that  as  total  flow 
through  a  capillary  bed  increases,  the  surface  area  and  volume  of  the  capillary 
bed  increase.  The  capillaries  in  any  region  are  fairly  uniformly  distributed  with 
their  flows  in  many  directions,  although  some  orientation  may  be  imposed  by 
the  shape  of  the  cells  of  a  region,  as  in  skeletal  muscle  where  most  of  the  capil¬ 
laries  run  parallel  to  the  large  striated  muscle  cells.  The  number  of  capillaries 
through  which  blood  is  flowing  may  vary  considerably,  depending  on  the  state 
of  the  tissue.  Capillary  flow  in  many  regions  is  intermittent  so  that  there  is  a 
rotation  of  function  among  the  various  capillaries  of  the  region.  Furthermore, 
there  are  many  interconnections  between  capillaries  so  that  one  cannot  think  in 
terms  of  many  single  channels  but  must  think  in  terms  of  a  meshwork  arising 
from  a  precapillary  arteriole  and  A-V  bridge  and  eventually  draining  into  a 
venule.  This  meshwork  is  further  interconnected  with  the  meshworks  arising 
from  other  precapillary  arterioles.  In  many  tissues  one  also  finds  arteriovenous 
shunts  which  bypass  the  capillary  beds. 

Although  there  is  some  variation  in  capillary  dimensions,  the  average  dimen¬ 
sions  commonly  quoted  in  the  literature  are  a  radius  of  .0004  cm.  and  a  length 
of  .04  cm.  The  blood  flowing  through  the  capillaries  consists  of  a  suspension  of 
cells  in  a  fluid,  the  plasma.  The  exchanges  we  will  be  concerned  with  will  be  pri¬ 
marily  between  the  plasma  and  extracellular  space,  so  we  will  speak  in  terms  of 
the  plasma  flow  to  a  region. 

To  complicate  our  problem  further,  we  must  consider  the  effects  of  recircula- 
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tion  and  of  mixing  in  the  heart  and  vessels.  Normally,  blood  flow  in  most  vessels 
involves  an  almost  uniform  velocity  flow  of  a  central  core  containing  the  sus¬ 
pended  cells.  The  major  part  of  the  shear  occurs  in  a  peripheral  zone  of  suspend¬ 
ing  medium,  the  plasma.  This  is  much  like  a  volume  displacement  flow.  We  shall 
therefore  assume  a  volume  displacement  flow  with  complete  mixing  in  a  direction 
perpendicular  to  the  flow  through  a  vessel.  Further  mixing  will  also  occur  in  the 
heart.  This  allows  us  to  concentrate  on  the  effects  introduced  by  the  time  lags 
inherent  in  the  recirculation. 

3.  A  mathematical  model:  general  considerations 

In  view  of  these  considerations,  what  sort  of  models  of  the  processes  occurring 
in  the  capillary  bed  might  we  consider? 

Right  at  the  outset  we  can  discard  the  possibility  of  setting  up  the  complete 
deterministic  equations  for  a  region.  This  would  require  complete  knowledge  of 
the  geometrical  relationships  between  all  cells,  extracellular  fluid  and  all  capil¬ 
laries  in  the  region.  Even  were  this  available  for  one  particular  area  of  one  par¬ 
ticular  animal  it  would  hardly  lead  to  tractable  equations.  One  might  consider 
using  average  values  for  cell  size,  intercellular  distance,  capillary  dimensions, 
and  intercapillary  distance,  but  this  would  not  make  the  problem  any  easier. 
Another  possibility,  which  has  been  suggested,  is  to  treat  a  region  as  consisting 
of  a  group  of  domains,  each  domain  consisting  of  a  cylindrical  region  surrounding 
a  capillary.  Aside  from  the  problem  of  choosing  appropriate  boundary  conditions 
and  solving  the  complex  partial  differential  equations  involved,  there  is  the  basic 
objection  that  this  may  be  an  inappropriate  model.  As  was  pointed  out  in  the 
previous  section,  the  capillary  bed  must  be  viewed  as  a  meshwork  of  intercon¬ 
necting  capillaries  and  not  a  series  of  single  channels  connecting  arterioles  and 
venules.  Since  point  by  point  measurements  of  the  concentration  of  a  compound 
in  extracellular  space  and  plasma  are  unavailable,  all  of  these  formulations  seem 
needlessly  complex.  A  macroscopic  viewpoint  may  be  more  useful.  For  polar 
molecules  we  assume  a  uniform  concentration  in  the  extracellular  space  and  in 
the  intracellular  space.  These  assumptions  are  justified  by  the  following  consid¬ 
erations:  at  least  from  a  macroscopic  viewpoint  the  fairly  uniform  distribution 
of  capillaries  and  capillary  flows  and  the  intermittency  of  capillary  action  tends 
to  smooth  out  diffusion  gradients  in  the  extracellular  space.  Furthermore,  the 
diffusion  of  polar  molecules  across  cell  membranes  into  the  intracellular  space  is 
much  slower  than  diffusion  across  the  capillary  membrane  into  the  extracellular 
space,  and  the  latter  is  much  slower  than  diffusion  within  the  extracellular  space 
[11],  [12],  [13].  In  essence  then,  we  use  a  spatial  average  for  the  extracellular 
concentration  and  the  above  considerations  at  least  give  us  some  assurance  that 
the  deviations  from  this  mean  value  will  not  oscillate  violently  as  one  moves 
from  point  to  point  in  the  extracellular  space. 

It  remains  then  to  discuss  how  one  can  treat  the  capillary  bed  in  view  of  the 
above  assumptions.  If  the  flow  rate  and  cross  sectional  area  for  the  various  capil- 
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laries  of  the  mesh  have  narrow  distributions,  we  can  lump  the  capillaries  of  a 
mesh  into  one  capillary  with  the  total  surface  area  and  plasma  volume  of  the 
mesh  and  a  length  determined  by  the  mean  transit  time  for  the  mesh. 

Thus  our  model  of  a  capillary  bed  consists  of  a  number  of  fairly  simple  func¬ 
tional  units,  each  unit  being  the  capillary  mesh  arising  from  a  precapillary 
arteriole.  These  units  are  connected  to  the  circulation  in  a  parallel  arrangement. 
Having  obtained  the  solution  to  the  equations  for  one  unit,  we  then  need  only 
sum  over  all  the  units  in  the  region  under  consideration  in  order  to  find  the 
amount  of  drug  in  intercellular  space,  extracellular  space  and  the  amount  of 
drug  entering  and  leaving  the  region  via  the  circulation.  For  the  time  being, 
however,  we  defer  consideration  of  this  last  statistical  aspect  of  the  problem, 
since  it  will  needlessly  complicate  matters  in  this  stage  of  the  investigation.  For 
the  time  being  we  commit  ourselves  to  the  use  of  a  mean  value  model. 


4.  Notation 

To  reduce  some  of  the  foregoing  ideas  to  mathematical  form,  let  us  introduce 
the  following  definitions : 

u(x,  t )  =  concentration  of  drug  in  moles  per  unit  volume  at  x  in  the  capillary 
at  time  t. 

v(t)  =  mean  concentration  of  drug  at  time  t  in  extracellular  space. 

w{t )  =  concentration  of  free  drug  in  intracellular  space. 

z{t)  =  concentration  of  drug-enzyme  complex  ED  resulting  from  reaction 

ki 

E  +  D^± ED  in  intracellular  space. 

ki 

E0  =  total  concentration  of  enzyme  in  intracellular  space. 
ki,  /c2  =  rate  constants  for  reaction  of  drug  with  E. 
ke  =  permeability  constant  for  capillary  walls. 
ki  =  permeability  constant  for  cell  membrane. 

K  =  transport  constant.  Net  active  transport  into  the  intracellular  space  is 
assumed  to  be  given  by  KA  lv.  This  is  only  an  approximation  for  most  active 
transport  systems. 

c  =  volume  flow  rate  of  plasma  in  capillary  bed. 

I  =  mean  capillary  length. 

Rp  =  plasma  volume  of  capillary  bed. 

Re  =  volume  of  extracellular  space. 

Ri  =  volume  of  intracellular  space. 


6.  A  model  of  a  capillary  bed 

Following  the  remarks  of  section  3,  we  may  diagram  our  model  of  the  capillary 
bed  as  shown  in  figure  1. 
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Figure  1 

Model  of  capillary  bed. 

u(x,  t)  —  concentration  of  drug  in  moles  per  unit  volume  at  x  in  the 
capillary  at  time  t, 

v(t )  =  mean  concentration  of  drug  at  time  t  in  extracellular  space, 
w(t)  =  concentration  of  free  drug  in  intracellular  space, 
z(0  =  concentration  of  drug-enzyme  complex  ED  resulting  from 
fei 

reaction  E  +  D  ^  ED  in  intracellular  space. 

Jh 


By  standard  arguments  we  obtain  the  set  of  equations 


(1) 


Hipj  ditj 


_ -  _  _r 

i  dt  3  dx 


du<  k.A( 


l 


(Uj  -  Vj), 


L 


Uj(x,  t )  dx  “  keAejVj  —  kiAij(Vj  —  Wj)  —  KjAijVj 


/r%\  n  dVj  keA  1 

<2)  K"rfT  =  “T 

(3)  Rij  ^  KjA^Vj  +  k{A  i](vj  tuy)  kiRi/Wj^Eoj  Zj)  ~ l-  Rijk%Zjf 


(4) 


^  =  kiWj(Eoj  -  Zj)  -  k2Zj. 


The  subscript  j  is  present  since  we  wish  to  speak  of  more  than  one  capillary  bed. 

This  is  a  particularly  interesting  set  of  equations,  consisting  of  a  partial  dif¬ 
ferential  equation,  an  integrodifferential  equation,  and  two  nonlinear,  first-order, 
ordinary  differential  equations.  We  should  point  out  that  from  the  viewpoint  of 
the  experimenter,  we  do  not  need  u{x,  t)  but  only  the  values  of  w(0,  t)  and  u(l,  t). 
These  are  the  quantities  which  are  measurable  in  practice.  It  should  further  be 
added  that  these  equations  are  incomplete  as  they  stand;  we  must  still  add 
equations  which  describe  the  effects  of  the  circulation  and  the  time  course  of 
drug  administration.  The  former  will  link  the  output  and  input  concentrations 
to  the  capillary  bed  via  the  time  lags  involved  in  recirculation  and  the  mixing 
in  the  circulation. 


6.  A  model  of  a  two-organ  being 

As  an  introduction,  let  us  consider  a  simplified  organism  (figure  2)  which  con¬ 
sists  of  two  organs  connected  in  parallel  to  the  circulation,  which  is  maintained 
by  a  simplified  heart.  Each  of  the  boxes  in  figure  2  represents  an  organ.  The 
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processes  occurring  in  each  organ  are  governed  by  equations  (1)  through  (4), 
with  j  taking  on  the  values  1  and  2  respectively.  To  these  we  must  add  the  equa¬ 
tions  of  the  circulation  and  of  drug  injection.  Let  us  assume  that  the  drug  is 
given  intravenously.  For  simplicity  we  will  treat  the  injection  as  though  it  were 


Diagram  of  a  simplified  organism  with  its  two  organs 
connected  in  parallel  to  its  circulation. 


given  directly  at  the  entrance  to  the  heart,  then  the  concentration  entering  the 
heart  will  be 


(5) 


0, 

UR(t)  = 

+  CjUijl,  t  -  r)  +  c2u2(l,  t  -  r)i 


t  <  0, 
OgKr, 

t  ^  r, 


where  r  is  the  circulation  time  lag  from  the  exit  of  the  capillary  beds  to  the  en¬ 
trance  to  the  heart,  ui(t)  describes  the  concentration  of  drug  due  to  the  injection, 
and  the  final  term  is  the  drug  concentration  resulting  from  recirculation.  We  are 
assuming  that  t  is  the  same  for  both  organs. 

The  heart  is  treated  as  a  mixing  chamber  of  volume 


(6) 


V*  = 


VB 


[log  (1  +  Vr/Vs)Y 

where  VE  is  the  ejection  volume  and  Vr  the  residual  volume  of  the  heart  [14], 
then  the  concentration  of  drug  ul  in  the  plasma  leaving  the  heart  is  given  by 

0,  t  <  0, 


(7) 


uL{t)  =  *|  c  _ 


—  e~ct/v *  J  UR(y)e?v/v *  dy,  t  >  0. 


Finally,  thp  concentration  of  drug  entering  the  capillary  beds  will  be 


rn 


n  <  *  << 
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It  is  apparent  now  that  these  conditions  make  our  set  of  differential  equations 
a  set  of  differential  difference  equations. 

7.  Computations 

The  equations  obtained  so  far  present  enormous  analytical  difficulties.  Al¬ 
though  the  steady  state  may  be  obtained  on  physical  grounds  it  is  the  time  course 
during  the  early  phases  of  the  process  which  we  will  need  if  we  wish  to  examine 
such  problems  as  the  relative  value  of  direct  intra-arteriol  injection  of  a  chemo¬ 
therapeutic  agent. 

We  can  circumvent  some  of  these  difficulties  if  we  can  simplify  the  ^-depend¬ 
ence  in  equation  (1).  A  simple  approximation  to  make  is  to  assume  that  the 
concentration  along  the  capillary  will  be  linear  in  x.  With  this  assumption  equa¬ 
tions  (1)  and  (2)  reduce  to  equations  (9)  and  (10). 

(9)  +  _  C,(M,0  -  -  .,) 

(10)  Rti  Aji  =  k.A.i  (M'°  J M  -  »,)  -  -  Wj)  -  KjAtjVj. 

In  these  equations  uj0  and  uji  represent  the  concentrations  in  the  plasma  enter¬ 
ing  and  leaving  capillary  bed  j. 

This  approximation  will  obviously  be  poor  for  processes  of  duration  less  than 
the  capillary  mesh  mean  transit  time;  this  is  particularly  true  for  those  periods 
when  a  concentration  wave  is  entering  or  leaving  the  capillary  bed.  If  we  con¬ 
sider  the  start  of  the  process  when  a  concentration  wave  is  entering  the  capillary 
bed,  we  have 

(11)  «ii(  0)  =  vj( 0)  =0,  t  =  0, 

and  equation  (8)  reduces  to  equation  (11), 

<i2)  ^  -  i  0  -  H*] »-«»  -  ^ 

Thus  if  duj0(0)/dt  —  0  then  duji(0)/dt  <  0  if  keAe}  >  2 cy,  giving  negative 
values  for  Uji(t)  for  some  short  interval  at  the  start  of  the  process.  On  the  other 
hand,  if  duj0(0)/dt  >  0  then  du,i(Q)/dt  can  still  be  negative  for  some  values  of 
keAej  <  2c j.  This  model  has  been  programmed  for  the  IBM  704  [15]  and,  as 
expected,  has  given  negative  concentration  waves  in  the  plasma  which  were 
propagated  for  some  time  because  of  the  recirculation.  Except  for  plasma  con¬ 
centration  waves  which  change  very  slowly  and  for  the  condition  keAe /  <  2c j, 
this  is  a  poor  approximation. 

To  avoid  the  difficulties  introduced  by  the  above  approximation,  the  capillary 
bed  can  be  replaced  by  n  mixing  chambers  in  series,  each  of  volume  Rpj/n  as 
shown  in  figure  3.  This  leads  to  the  set  of  equations  (13)  and  equation  (14)  in 
place  of  equations  (1)  and  (2)  respectively. 
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/<  dujm  nCj  ,  .  keA ej  ,  >,  i  o 

(13)  ,  —  p  i  Ujm)  „  (ujm  Vj),  in  —  1,  2,  •  •  •  ,  n} 

U/l  pj  ii/pj 

(14)  R„  ^  f  (u,«  -  »,)  -  Mofe  -  to,)  -  KtAuv,. 

This  approximating  model  has  the  advantage  that  negative  solutions  cannot 
be  obtained  and  that  equations  (13)  and  (14)  converge  on  (1)  and  (2)  respec¬ 
tively  as  n  increases.  Furthermore,  we  expect  on  physical  grounds  that  we  need 
only  make  n  large  enough  such  that  RPj/n  is  less  than  the  volume  of  the  plasma 
occupied  by  the  concentration  wave  entering  the  capillary  bed,  to  obtain  a  good 


i  2  n 


Figure  3 

Diagram  of  n  mixing  chambers  in  series. 


approximation  to  equations  (1)  and  (2).  This  model  has  also  been  programmed 
for  the  IBM  704  [15]  and  has  given  reasonable  results  for  the  time  course  of 
concentration  changes  in  the  capillary  plasma,  extracellular  space,  and  intra¬ 
cellular  space.  However,  mixing  in  the  circulation  is  slower  than  would  be  ex¬ 
pected  from  experimental  data  on  equilibration  of  dyes  and  713,-labeled  albumin 
after  injection  into  the  circulation  [16],  [17].  This  is  probably  due  both  to  the 
assumption  of  volume  displacement  flow  and  the  use  of  only  two  organs  with 
equal  circulatory  time  lags.  We  would  expect  more  rapid  mixing  in  a  many-organ 
model  with  a  distribution  of  time  lags.  Although  the  required  mixing  in  the  cir¬ 
culation  could  be  introduced  arbitrarily  by  assuming  a  distribution  of  mixing 
chambers  in  the  circulation,  it  would  be  more  realistic  and  esthetically  more 
satisfying  to  incorporate  more  of  the  known  features  of  blood  flow  in  the  model 
of  the  circulation. 

8.  Discussion 

The  equations  obtained  from  the  model  of  drug  distribution  in  a  two-organ 
being  are  quite  complex  analytically.  We  would  expect  these  difficulties  to  be 
multiplied  manyfold  in  a  model  of  a  many-organ  being  which  would  be  required 
if  we  wished  to  approximate  the  situation  in  a  mammalian  organism.  Fortu¬ 
nately,  the  use  of  high-speed  computers  allows  us  to  circumvent  many  of  the 
analytical  difficulties.  Once  we  recognize  this  we  can  treat  many  biological 
problems  in  a  more  realistic  fashion. 

A  projection  of  these  concepts,  which  we  believe  will  be  of  major  importance 
for  biological  research,  is  the  idea  of  experimentation  with  the  computer.  How¬ 
ever,  this  is  a  useful  concept  only  after  fairly  realistic  models  of  a  biological 
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system  can  be  formulated.  Once  this  stage  has  been  reached  we  can  consider 
running  experiments  on  the  model  system  with  the  computer.  If  this  is  done  in 
conjunction  with  experiments  on  the  real  biological  system  one  obtains  a  con¬ 
tinuing  feedback  process  between  experiments  on  the  model  and  on  the  real 
system  which  can  lead  to  a  great  saving  in  time  and  experimental  effort.  Perhaps 
the  most  significant  gain  to  be  expected  from  such  an  approach  is  the  clarification 
of  basic  concepts  which  must  ensue  since  the  setting  up  of  the  model  system 
requires  a  rigorous  formulation  of  our  ideas  about  the  basic  mechanisms  at  work 
in  the  real  system  under  study.  Applying  these  ideas  to  a  model  of  drug  dis¬ 
tribution,  let  us  consider  that  we  have  a  class  of  drugs  whose  biological  effect  is 
primarily  due  to  a  reaction  with  one  component  of  the  intracellular  space.  By 
a  variation  of  parameters,  the  computer  can  then  be  used  to  study  the  effects 
of  blood  flow,  permeability  constants,  and  reaction  rate  constants  on  drug  action. 
This  would  then  provide  criteria  for  picking  or  redesigning  the  drug  for  certain 
optimal  effects.  Similarly,  experiments  could  be  run  to  design  optimal  policies 
of  drug  administration. 
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